59. 协方差分析

from pathlib import Path  # 导入Path以兼容两种课程运行平台
import numpy as np  # 导入NumPy用于矩阵与有限值计算
import pandas as pd  # 导入Pandas用于行情整理与关系计算
linux_root = Path('/home/ubuntu/r2_data_mount/data')  # 指定Linux课程数据根目录
windows_root = Path(r'C:\qiufei\data')  # 指定Windows课程数据根目录
data_root = linux_root if linux_root.exists() else windows_root  # 选择当前可用根目录
price_path = data_root / 'stock' / 'stock_price_pre_adjusted.h5'  # 定位前复权行情快照
basic_path = data_root / 'stock' / 'stock_basic_data.h5'  # 定位公司基础信息文件
target_codes = ['002415.SZ','600276.SH','600104.SH','002230.SZ']  # 定义四只长三角非金融公司样本
expected_names = {'002415.SZ':'海康威视','600276.SH':'恒瑞医药','600104.SH':'上汽集团','002230.SZ':'科大讯飞'}  # 定义待核验代码名称映射
is_price_available = price_path.exists()  # 记录行情数据要求状态
if not is_price_available:  # 行情文件缺失时提示读者下载
    print(f'未找到课程数据文件：{price_path}；请从课程数据下载入口获取后重新运行。')  # 提示读者核对文件位置
    print('数据读取成功后再计算滚动相关和组合方差。')  # 说明下一步

name_audit = {'status':'NOT_VERIFIED','observed_names':{}}  # 初始化名称未核验状态
if is_price_available and basic_path.exists():  # 公司基础信息存在时检查名称
    companies = pd.read_hdf(basic_path,key='stock_basic_info').rename(columns={'order_book_id':'ts_code','symbol':'name'})  # 读取并整理公司基础信息
    companies['ts_code'] = companies['ts_code'].str.replace('.XSHG','.SH',regex=False).str.replace('.XSHE','.SZ',regex=False)  # 统一证券代码后缀
    name_column = next((column for column in ['name','股票简称','sec_name'] if column in companies.columns),None)  # 选择可用名称字段
    if name_column is None or 'ts_code' not in companies.columns:  # 核验名称与代码字段要求
        raise ValueError('公司基本信息中缺少 ts_code 或公司名称字段')  # 提示缺少的字段
    observed_names = companies.loc[companies['ts_code'].isin(target_codes),['ts_code',name_column]].drop_duplicates('ts_code').set_index('ts_code')[name_column].to_dict()  # 提取目标代码名称
    mismatched_names = {code:observed_names.get(code) for code,name in expected_names.items() if observed_names.get(code)!=name}  # 识别不一致映射
    if mismatched_names:  # 名称不一致时停止使用真实公司标签
        raise ValueError(f'股票代码与公司名称不匹配：{mismatched_names}')  # 提示不一致的记录
    name_audit = {'status':'PASS','observed_names':observed_names}  # 保存通过的名称核验结果
elif is_price_available:  # 仅行情可用时保留未核验状态
    print(f'未找到公司基本信息文件：{basic_path}；请下载后重新核对公司名称。')  # 提示读者下载文件

if is_price_available:  # 行情存在时读取同一时间段真实价格
    rq_codes = [code.replace('.SH','.XSHG').replace('.SZ','.XSHE') for code in target_codes]  # 转换为HDF证券代码
    prices = pd.read_hdf(price_path,key='data',where='order_book_id in rq_codes',columns=['close']).reset_index().rename(columns={'order_book_id':'ts_code','date':'trade_date'})  # 读取目标前复权收盘价
    prices['ts_code'] = prices['ts_code'].str.replace('.XSHG','.SH',regex=False).str.replace('.XSHE','.SZ',regex=False)  # 恢复课堂证券代码后缀
    prices['trade_date'] = pd.to_datetime(prices['trade_date'])  # 统一交易日期类型
    sample = prices.loc[prices['ts_code'].isin(target_codes)&prices['trade_date'].between('2019-01-01','2025-12-31')]  # 固定样本资产与期间
    observed_codes = set(sample['ts_code'].dropna().unique())  # 收集实际到齐资产代码
    missing_codes = sorted(set(target_codes)-observed_codes)  # 识别缺失目标资产
    if missing_codes:  # 资产未全部到齐时停止矩阵计算
        raise ValueError(f'价格数据中缺少这些股票：{missing_codes}')  # 提示缺少的股票
    if sample.duplicated(['trade_date','ts_code']).any():  # 检查资产日期键唯一性
        raise ValueError('同一股票与交易日存在重复记录')  # 提示重复记录
    raw_close_panel = sample.pivot(index='trade_date',columns='ts_code',values='close').reindex(columns=target_codes)  # 构造固定顺序价格面板
    missing_observations = raw_close_panel.isna().sum()  # 统计各资产缺失观测
    close_panel = raw_close_panel.dropna()  # 仅保留共同交易日
    if close_panel.empty:  # 无共同交易日时停止分析
        raise ValueError('这些股票之间没有共同交易日，请核对日期范围')  # 提示核对日期

if is_price_available:  # 数据就绪时计算收益与关系矩阵
    return_panel = close_panel.pct_change().dropna()  # 计算共同交易日日收益率
    covariance_matrix = return_panel.cov()  # 计算全样本协方差矩阵
    correlation_matrix = return_panel.corr()  # 计算全样本相关矩阵
    if len(return_panel)<80:  # 保证60日窗口与四周满足条件所需长度
        raise ValueError('共同交易日不足，无法评估60日窗口与连续四周满足条件')  # 样本不足时停止

def summarize_weekly_pair(returns_frame,left_code,right_code):  # 定义单资产对周度状态函数
    rolling_pair = returns_frame[left_code].rolling(60).corr(returns_frame[right_code]).dropna()  # 计算60日滚动相关
    weekly_pair = rolling_pair.resample('W-FRI').last().dropna()  # 转为周末最后状态
    if len(weekly_pair)<4:  # 核验连续四周所需样本
        raise ValueError(f'{left_code}/{right_code}不足四周数据')  # 提示样本不足
    weekly_breach = weekly_pair.gt(0.8)  # 标记正相关超课堂阈值状态
    breach_streak = weekly_breach.groupby((~weekly_breach).cumsum()).cumsum().astype(int)  # 计算连续满足条件周数
    pair_audit = pd.DataFrame({'asset_pair':f'{left_code} / {right_code}','week_end':weekly_pair.index,'correlation_60d':weekly_pair.values,'positive_concentration':weekly_breach.values,'continuous_weeks':breach_streak.values})  # 生成完整周度检查表
    current_state = {'asset_pair':f'{left_code} / {right_code}','week_end':weekly_pair.index[-1],'correlation_60d':weekly_pair.iloc[-1],'positive_threshold':.8,'continuous_weeks':int(breach_streak.iloc[-1]),'four_week_condition_met':bool(breach_streak.iloc[-1]>=4)}  # 生成当前资产对高正相关状态
    return pair_audit.tail(4).to_dict('records'),current_state  # 返回最近四周检查与当前状态

if is_price_available:  # 数据就绪时汇总全部资产对
    weekly_rows = []  # 初始化当前状态容器
    weekly_audit_rows = []  # 初始化最近四周检查容器
    for left_position,left_code in enumerate(return_panel.columns):  # 逐个选择左侧资产
        for right_code in return_panel.columns[left_position+1:]:  # 仅遍历不重复右侧资产
            pair_audit_rows,pair_state = summarize_weekly_pair(return_panel,left_code,right_code)  # 计算单资产对状态
            weekly_audit_rows.extend(pair_audit_rows)  # 累计资产对最近四周检查
            weekly_rows.append(pair_state)  # 累计资产对当前状态
    weekly_state_table = pd.DataFrame(weekly_rows)  # 汇总当前周度状态表
    weekly_audit_table = pd.DataFrame(weekly_audit_rows)  # 汇总最近四周检查表
    if len(weekly_state_table)!=6 or weekly_state_table['asset_pair'].nunique()!=6:  # 核验四资产全部六个唯一资产对
        raise ValueError('输入文件、字段、样本量或数值不符合当前分析要求，请按本页说明检查')  # 资产对覆盖不足时停止

if is_price_available:  # 数据就绪时计算课堂组合风险
    portfolio_weights = pd.Series([.25,.25,.25,.25],index=target_codes,name='weight').reindex(return_panel.columns)  # 建立等权课堂组合并对齐资产顺序
    if len(portfolio_weights)!=return_panel.shape[1] or portfolio_weights.isna().any():  # 核验权重长度与资产顺序
        raise ValueError('权重数量或股票顺序与收益率矩阵不一致')  # 提示核对权重
    if (portfolio_weights<0).any() or not np.isclose(portfolio_weights.sum(),1.0):  # 核验非负与和为1
        raise ValueError('权重必须为非负数且合计为 1')  # 提示权重要求
    annual_covariance = return_panel.tail(60).cov()*252  # 用课堂252因子年化最近60日协方差
    portfolio_variance = float(portfolio_weights@annual_covariance@portfolio_weights)  # 计算组合年化方差
    risk_budget = .04  # 设置年化20%波动对应的课堂方差预算假设
    correlation_condition_met = bool(weekly_state_table['four_week_condition_met'].any())  # 判断是否存在连续四周高正相关
    variance_condition_met = bool(portfolio_variance>risk_budget)  # 判断组合方差是否超过课堂预算
    portfolio_condition_met = correlation_condition_met or variance_condition_met  # 合并两类核对判断条件

if is_price_available:  # 数据就绪时整理结果并比较不同设置
    high_correlation_pairs = weekly_state_table.loc[weekly_state_table['four_week_condition_met'],'asset_pair'].tolist()  # 提取当前连续四周高正相关资产对
    decision_table = pd.DataFrame([{'state_window_days':60,'state_sampling':'W-FRI','sensitivity_windows_weeks':'20/60/120','annualization_factor':252,'variance_unit':'年化日收益方差','portfolio_variance':portfolio_variance,'risk_budget':risk_budget,'positive_threshold':.8,'correlation_condition':'correlation_60d > 0.8 连续四周','high_correlation_pairs':high_correlation_pairs,'correlation_condition_met':correlation_condition_met,'variance_condition_met':variance_condition_met,'condition_met':portfolio_condition_met,'owner':'风控负责人','action':'核对共同因子并降低集中权重' if portfolio_condition_met else '保持权重并按周观察','review_note':'全部资产对不再连续四周高于0.8且年化方差回到0.04内；下一季度再次检查'}])  # 整理课堂结果和建议
    weekly_return_panel = close_panel.resample('W-FRI').last().pct_change(fill_method=None).dropna()  # 从真实周末价格生成周收益，避免把交易日误称为周
    sensitivity_rows = []  # 初始化周窗口与阈值敏感性结果
    for window_weeks in [20,60,120]:  # 比较三种事先设定周收益窗口
        for positive_threshold in [.7,.8,.9]:  # 比较三档课堂正相关阈值
            qualifying_pair_count = 0  # 统计曾满足连续四周条件的资产对数
            qualifying_week_count = 0  # 统计全部资产对的满足条件周数
            for left_position,left_code in enumerate(weekly_return_panel.columns):  # 逐个选择左侧资产
                for right_code in weekly_return_panel.columns[left_position+1:]:  # 仅遍历不重复右侧资产
                    rolling_weekly_correlation = weekly_return_panel[left_code].rolling(window_weeks).corr(weekly_return_panel[right_code])  # 按周收益窗口计算相关序列
                    weekly_threshold_breach = rolling_weekly_correlation.gt(positive_threshold)  # 标记正相关越过课堂阈值的周
                    four_week_condition_met = weekly_threshold_breach.rolling(4,min_periods=4).sum().eq(4)  # 仅连续四周越界时满足条件
                    qualifying_pair_count += int(four_week_condition_met.any())  # 汇总曾满足条件的资产对
                    qualifying_week_count += int(four_week_condition_met.sum())  # 汇总持续满足条件的周数
            sensitivity_rows.append({'window_weeks':window_weeks,'positive_threshold':positive_threshold,'pair_count':6,'qualifying_pairs':qualifying_pair_count,'qualifying_weeks':qualifying_week_count})  # 保存覆盖全部6对的频率清晰敏感性情景
    sensitivity_table = pd.DataFrame(sensitivity_rows)  # 汇总九种周窗口与阈值情景
    if len(sensitivity_table)!=9 or not sensitivity_table['pair_count'].eq(6).all():  # 核验九情景均覆盖四资产全部6对
        raise ValueError('输入文件、字段、样本量或数值不符合当前分析要求，请按本页说明检查')  # 敏感性覆盖不足时停止
    student_portfolio_variance,student_weekly_state_table,student_sensitivity_table = portfolio_variance,weekly_state_table.copy(),sensitivity_table.copy()  # 同时锁定本章计算方差、全部6对状态与九行情景以免后续示例数据覆盖
    strongest_pair = correlation_matrix.where(~np.eye(len(correlation_matrix),dtype=bool)).stack().idxmax()  # 定位最高样本相关资产对
    print({'status':'READY','name_audit':name_audit['status'],'assets':len(observed_codes),'rows':len(return_panel),'period':[str(return_panel.index.min().date()),str(return_panel.index.max().date())],'missing':missing_observations.to_dict()})  # 输出绑定期间与覆盖依据
    print({'strongest_pair':' / '.join(strongest_pair),'correlation':round(float(correlation_matrix.loc[strongest_pair]),3),'portfolio_variance':round(portfolio_variance,6),'risk_budget':risk_budget,'condition_met':portfolio_condition_met})  # 输出关系与满足条件依据
    print({'weekly_pairs':len(weekly_state_table),'sensitivity_scenarios':len(sensitivity_table),'owner':decision_table.at[0,'owner'],'action':decision_table.at[0,'action']})  # 输出比较范围和建议

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['Source Han Serif SC']
plt.rcParams['axes.unicode_minus'] = False

# 创建中性多资产收益率纯教学示例，不承载企业或投资依据
np.random.seed(42)
n_days = 500

# 各资产的期望日收益率
mean_returns = [0.001, 0.0008, 0.0012, 0.0005]

# 设定教学示例协方差矩阵，非市场观测
cov_matrix_true = np.array([
    [0.0004, 0.0002, 0.00015, 0.0001],
    [0.0002, 0.0003, 0.0001, 0.00005],
    [0.00015, 0.0001, 0.00035, 0.00008],
    [0.0001, 0.00005, 0.00008, 0.0002]
])

# 生成多元正态分布数据
returns = np.random.multivariate_normal(mean_returns, cov_matrix_true, n_days)

# 转换为DataFrame
assets = ['资产A', '资产B', '资产C', '资产D']  # 使用中性资产标签避免企业依据误读
df_returns = pd.DataFrame(returns, columns=assets)

print('教学示例，非市场观测｜收益率数据(前10行):')  # 在输出表前声明教学示例边界
print(df_returns.head(10))

# 计算协方差矩阵
cov_matrix = df_returns.cov()
print('教学示例，非市场观测｜协方差矩阵:')  # 在人工协方差输出前声明依据边界
print(cov_matrix.round(6))

# 计算Pearson相关系数矩阵
corr_matrix = df_returns.corr()
print('\n教学示例，非市场观测｜相关系数矩阵:')  # 在人工相关矩阵输出前声明依据边界
print(corr_matrix.round(4))

# 提取上三角矩阵（排除对角线），避免重复
mask = np.triu(np.ones_like(corr_matrix, dtype=bool), k=1)
corr_upper = corr_matrix.where(mask)

# 找最大和最小相关系数
max_corr = corr_upper.max().max()
min_corr = corr_upper.min().min()
max_pair = corr_upper.stack().idxmax()
min_pair = corr_upper.stack().idxmin()

print(f'教学示例，非市场观测｜最高相关性: {max_pair} = {max_corr:.4f}')  # 在孤立输出保留依据边界
print(f'最低相关性: {min_pair} = {min_corr:.4f}')

max_sharpe_idx = (df_portfolios['收益率'] / df_portfolios['风险(标准差)']).idxmax()
optimal_weights = weights_list[max_sharpe_idx]

print('教学示例，非市场观测｜样本内最高收益风险比权重:')  # 输出示例数据内描述结果而非投资建议
for asset, weight in zip(assets_list, optimal_weights):
    print(f'{asset}: {weight:.2%}')

print('教学示例，非市场观测｜相关系数统计:')  # 在孤立输出表头保留依据边界
print(f'均值: {rolling_corr.mean():.4f}')
print(f'标准差: {rolling_corr.std():.4f}')
print(f'最大值: {rolling_corr.max():.4f}')
print(f'最小值: {rolling_corr.min():.4f}')

# 识别相关系数跌破课堂阈值的时期；未做统计检验，不称显著
threshold = rolling_corr.mean() - 2 * rolling_corr.std()
low_corr_periods = rolling_corr[rolling_corr < threshold]
print(f'\n相关性异常低的时期({len(low_corr_periods)}天):')
print(low_corr_periods.head())

high_corr_pairs = []
for i in range(n_industries):
    for j in range(i+1, n_industries):
        corr_val = industry_corr.iloc[i, j]
        if abs(corr_val) > 0.5:
            high_corr_pairs.append((industries_network[i], industries_network[j], corr_val))

print('教学示例，非市场观测｜超过课堂阈值的中性类别对(|ρ| > 0.5):')  # 输出带依据边界的阈值检索结果
for pair in sorted(high_corr_pairs, key=lambda x: abs(x[2]), reverse=True):
    print(f'{pair[0]} - {pair[1]}: {pair[2]:.2f}')

from scipy.stats import pearsonr
from sklearn.linear_model import LinearRegression

np.random.seed(42)
n = 500

# 构造三变量数据
market = np.random.normal(0.001, 0.02, n)
industry_factor = np.random.normal(0, 0.01, n)
stock = 0.8 * market + 0.5 * industry_factor + np.random.normal(0, 0.015, n)

df_partial = pd.DataFrame({
    '市场': market,
    '个股': stock,
    '行业因子': industry_factor
})

# 简单相关系数
corr_market_stock = df_partial['市场'].corr(df_partial['个股'])
corr_industry_stock = df_partial['行业因子'].corr(df_partial['个股'])

print('教学示例，非市场观测｜简单相关系数:')  # 在孤立输出表头保留依据边界
print(f'市场-个股: {corr_market_stock:.4f}')
print(f'行业因子-个股: {corr_industry_stock:.4f}')

# 个股对行业因子回归，取残差
model_stock = LinearRegression().fit(
    df_partial[['行业因子']], df_partial['个股']
)
residual_stock = df_partial['个股'] - model_stock.predict(df_partial[['行业因子']])

# 市场对行业因子回归，取残差
model_market = LinearRegression().fit(
    df_partial[['行业因子']], df_partial['市场']
)
residual_market = df_partial['市场'] - model_market.predict(df_partial[['行业因子']])

# 偏相关系数 = 残差的相关系数
partial_corr = np.corrcoef(residual_market, residual_stock)[0, 1]

print(f'教学示例，非市场观测｜偏相关系数(控制行业因子):')  # 在孤立输出表头保留依据边界
print(f'市场-个股: {partial_corr:.4f}')
print(f'解释: 剔除示例数据中第三变量的线性部分后，计算两组残差的样本相关性')  # 输出教学示例的统计边界而非业务结论